p8105_final_project

Introduction and Motivation

The five leading causes of death in the United States from 1999 to 2014 are cancer, heart disease, unintentional injury, chronic lower respiratory disease, and stroke. The dataset includes the U.S. Department of Health and Human Services public health regions. Therefore, we can investigate the leading causes of death of each region, and develop accordingly public health policy and remedies.

Data Descriptions

• Version of data used for analysis is Here.
  

Approaches

• From the dataset, we obtain the information that the number of potentially excess deaths from the five leading causes in rural areas was higher than those in urban areas.

• We then analyzed several factors that might influence the rural-urban difference in potentially excess deaths from the five leading causes, many of which are associated with sociodemographic and ecological differences between rural and urban areas.

• Through statistical analysis, our report provides an interactive and straightforward view on the potentially excess deaths from the five leading causes of death in non-metropolitan and metropolitan areas.

• The ultimate goal is to bring attention to preventing deaths in the rural areas through improving healthcare services and public health programs.

Visualizations

• Link to the Dashboard for this project is Here.
• Link to the Shiny Dashboard for this project is Here.

Conclusion

Notes

• Link to the GitHub repository containing R Markdown document for this project is Here.
• Link to the Screencast for this project is Here.
• Link to the GitHub repository containing R Markdown document for Shiny Dashboard is Here.

cod_data = read_csv("./data/NCHS_-_Potentially_Excess_Deaths_from_the_Five_Leading_Causes_of_Death.csv") %>%
  clean_names() %>%
  na.omit() %>%
  filter(!(state == "United States")) %>%
  separate(., percent_potentially_excess_deaths, into = c("percent_excess_death"), sep = "%") %>% 
  mutate(percent_excess_death = as.numeric(percent_excess_death), mortality = observed_deaths/population * 10000, mortality = as.numeric(mortality)) %>% 
  select(year, age_range, cause_of_death, state, locality, observed_deaths, population, expected_deaths, potentially_excess_deaths, percent_excess_death, mortality, hhs_region)

## Parsed with column specification:
## cols(
##   Year = col_integer(),
##   `Cause of Death` = col_character(),
##   State = col_character(),
##   `State FIPS Code` = col_character(),
##   `HHS Region` = col_integer(),
##   `Age Range` = col_character(),
##   Benchmark = col_character(),
##   Locality = col_character(),
##   `Observed Deaths` = col_integer(),
##   Population = col_integer(),
##   `Expected Deaths` = col_integer(),
##   `Potentially Excess Deaths` = col_integer(),
##   `Percent Potentially Excess Deaths` = col_character()
## )

## Warning: Too many values at 191748 locations: 1, 2, 3, 4, 5, 6, 7, 8, 9,
## 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...

##columns removed

 #"state_fips_code"                      "benchmark"   "potentially_excess_deaths" "percent_excess_death"      "mortality"   

• the year variable contains data collected from 2005-2015.
• filtered out United States in the state variable.
• Added a variable mortality which is calculated by observed_deaths/population * 10000. This variable indicates the number of deathes observed in every 10000 people in the three geographic regions: Metropolitan, Nonmetropolitan and All.

xs2329:

cod_data %>%
  group_by(cause_of_death) %>% 
  ggplot(aes(x = locality, y = mortality, fill = cause_of_death)) +
  geom_col() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title="Locality vs. Mortality") +
  labs(x="locality", y="mortality") 

This bar graph shows the distribution of cause of death within mortality in the three geographic regions: Metropolitan, Nonmetropolitan and All. We observe that the weight of each composition of cuases of death is the same regardless of the locality. Nonmetropolitan area seems to have the highest mortaility, the number of deathes observed in every 10000 people.

cod_data %>%
  group_by(cause_of_death) %>% 
  ggplot(aes(x = locality, y = mortality)) + geom_boxplot(aes(color = cause_of_death), na.rm = T) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title="Locality vs. Mortality") +
  labs(x="locality", y="mortality") 

We observe that the median of mortaility is the highest among all 5 causes of death, followed by heart disease, unintentional injury, heart disease, chronic lowee respiratory disease and stroke. This pattern remains consistant when we subdivide locaility into metropolitan region and nonmetrpolitan region, which is identical to our previous finding. All 5 cuases of death are right skewed regardless of the locality.

mortality_lm = lm(mortality ~locality, data = cod_data)
summary(mortality_lm)

## 
## Call:
## lm(formula = mortality ~ locality, data = cod_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1454 -2.9547 -1.2983  0.9988 19.8959 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              4.10331    0.01672 245.474   <2e-16 ***
## localityMetropolitan    -0.22576    0.02372  -9.519   <2e-16 ***
## localityNonmetropolitan  1.14451    0.02438  46.942   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.317 on 191745 degrees of freedom
## Multiple R-squared:  0.01821,    Adjusted R-squared:  0.0182 
## F-statistic:  1778 on 2 and 191745 DF,  p-value: < 2.2e-16

interpretation on linear regression:

Here we obtain the linear regression model of mortality = 4.103 - 0.226localityMetropolitan -1.145localityNonmetropolitan. We expect to see 0.22576 less deaths for every 10000 people in the metropolitan region as compared to all regions. We expect to see 1.14451 more death for every 10000 people in the nonmetropolitan region as compared to all regions.

cod_data %>%
  mutate(year = as.factor(year))

## # A tibble: 191,748 x 12
##      year age_range cause_of_death   state        locality observed_deaths
##    <fctr>     <chr>          <chr>   <chr>           <chr>           <int>
##  1   2005      0-49         Cancer Alabama             All             756
##  2   2005      0-49         Cancer Alabama    Metropolitan             556
##  3   2005      0-49         Cancer Alabama Nonmetropolitan             200
##  4   2005      0-49         Cancer Alabama             All             756
##  5   2005      0-49         Cancer Alabama    Metropolitan             556
##  6   2005      0-49         Cancer Alabama Nonmetropolitan             200
##  7   2005      0-49         Cancer Alabama             All             756
##  8   2005      0-49         Cancer Alabama    Metropolitan             556
##  9   2005      0-49         Cancer Alabama Nonmetropolitan             200
## 10   2005      0-54         Cancer Alabama             All            1346
## # ... with 191,738 more rows, and 6 more variables: population <int>,
## #   expected_deaths <int>, potentially_excess_deaths <int>,
## #   percent_excess_death <dbl>, mortality <dbl>, hhs_region <int>

p <- cod_data %>%
  plot_ly(
    x = ~expected_deaths, 
    y = ~observed_deaths, 
    size = ~population, 
    color = ~cause_of_death, 
    frame = ~hhs_region, 
    text = ~state, 
    hoverinfo = "text",
    type = 'scatter',
    mode = 'markers'
  ) 
p

### Question of concern: If percent excess death has significant difference between metro and non-metro groups.
gp_cod_data = cod_data %>%
  group_by(year, locality) %>% 
  summarise(mean_ped = mean(percent_excess_death))

ggplot(gp_cod_data, aes(x = year, y = mean_ped, fill = locality)) + geom_point(stat = "identity") + geom_smooth() +
facet_grid(. ~ locality) 

## `geom_smooth()` using method = 'loess'

gp_locality_lm = lm(mean_ped ~locality , data = gp_cod_data)
summary(gp_locality_lm)

## 
## Call:
## lm(formula = mean_ped ~ locality, data = gp_cod_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5667 -1.0325 -0.3470  0.7922  3.2937 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              33.8866     0.4071  83.238  < 2e-16 ***
## localityMetropolitan     -2.0599     0.5757  -3.578   0.0012 ** 
## localityNonmetropolitan   8.0460     0.5757  13.975 1.13e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.35 on 30 degrees of freedom
## Multiple R-squared:  0.9198, Adjusted R-squared:  0.9145 
## F-statistic: 172.1 on 2 and 30 DF,  p-value: < 2.2e-16

broom::tidy(gp_locality_lm)

##                      term  estimate std.error statistic      p.value
## 1             (Intercept) 33.886633 0.4071070 83.237648 4.783207e-37
## 2    localityMetropolitan -2.059931 0.5757363 -3.577906 1.200041e-03
## 3 localityNonmetropolitan  8.045979 0.5757363 13.975111 1.132457e-14

The regression model is mean percent_excess_death = 33.89 - 2.05*Metropolitan + 8.046Nonmetropolitan. For people in metropolitan, the expected percent excess death has a decrease of 2.06. For people in non-metropolitan aream the expected percent excess death has an increase of 8.046. The p-values show these predictors are both significant. The adjusted R squared is 0.91, which indicates the regression model explains 91% of variation in mean percent excess death due to the variation in locality.

# How should we help nonmetropolitan areas improve public health? 
cod_datayr2 = read_csv("./data/NCHS_-_Potentially_Excess_Deaths_from_the_Five_Leading_Causes_of_Death.csv") %>%
  clean_names() %>%
  na.omit() %>%
  filter(!(state == "United States")) %>%
  separate(., percent_potentially_excess_deaths, into = c("percent_excess_death"), sep = "%") %>%
  mutate(percent_excess_death = as.numeric(percent_excess_death), mortality = observed_deaths/population * 10000, mortality = as.numeric(mortality)) %>%
  select(year, cause_of_death, state, benchmark, locality, observed_deaths, population, expected_deaths, potentially_excess_deaths, percent_excess_death, mortality)

## Parsed with column specification:
## cols(
##   Year = col_integer(),
##   `Cause of Death` = col_character(),
##   State = col_character(),
##   `State FIPS Code` = col_character(),
##   `HHS Region` = col_integer(),
##   `Age Range` = col_character(),
##   Benchmark = col_character(),
##   Locality = col_character(),
##   `Observed Deaths` = col_integer(),
##   Population = col_integer(),
##   `Expected Deaths` = col_integer(),
##   `Potentially Excess Deaths` = col_integer(),
##   `Percent Potentially Excess Deaths` = col_character()
## )

## Warning: Too many values at 191748 locations: 1, 2, 3, 4, 5, 6, 7, 8, 9,
## 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...

average percent access death for each cause in different locality 2005-2015

cod_datayr2 %>%
  filter(benchmark == "2010 Fixed")%>%
  group_by(cause_of_death,year, locality)%>%
  mutate(year_percent_death=mean(percent_excess_death))%>%
  distinct(year, cause_of_death, locality, .keep_all = TRUE)%>%
  ungroup(year, locality)%>%
ggplot(aes(x = year, y = year_percent_death, color = cause_of_death))+
  geom_point(size = 1)+
  geom_line(size = 1) +
  facet_grid(.~locality)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  scale_x_continuous( breaks=seq(2005,2015,1))+
  labs(title = "Average yearly excess death percentage for different locality")

Interpretation: As can be seen in the plot, if we stratify by cause, we will find some increasing pattern in subgroups(unintentional injury) in certain period of time(2010-2015). In nonmetropolitan regions, unintentional injury and chronic lower respiratory disease have on average the highest excess death percentage in nonmetropolitan areas and needs more attention from public health experts.

cod_datayrlm = cod_datayr2 %>%
  filter(locality == "Nonmetropolitan",
         benchmark == "2010 Fixed")
lmpop = lm(data=cod_datayrlm, percent_excess_death~year+cause_of_death)
summary(lmpop)

## 
## Call:
## lm(formula = percent_excess_death ~ year + cause_of_death, data = cod_datayrlm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -55.456  -9.828   0.576  10.404  39.489 
## 
## Coefficients:
##                                                  Estimate Std. Error
## (Intercept)                                     769.06335   67.48938
## year                                             -0.36941    0.03358
## cause_of_deathChronic Lower Respiratory Disease  27.06490    0.33901
## cause_of_deathHeart Disease                      13.90257    0.32844
## cause_of_deathStroke                             12.91441    0.33784
## cause_of_deathUnintentional Injury               31.18470    0.32812
##                                                 t value Pr(>|t|)    
## (Intercept)                                       11.39   <2e-16 ***
## year                                             -11.00   <2e-16 ***
## cause_of_deathChronic Lower Respiratory Disease   79.83   <2e-16 ***
## cause_of_deathHeart Disease                       42.33   <2e-16 ***
## cause_of_deathStroke                              38.23   <2e-16 ***
## cause_of_deathUnintentional Injury                95.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.92 on 19718 degrees of freedom
## Multiple R-squared:  0.365,  Adjusted R-squared:  0.3649 
## F-statistic:  2267 on 5 and 19718 DF,  p-value: < 2.2e-16

From the linear model based on 2010 fixed benchmark, we can tell that controlling for cause of death, the excess death percentage is expected to decrease by 0.37 each year. Overall, the excess death percentage shows a decreasing pattern with time. Besides, in the same year, compared with cancer, other 4 causes have significantly higher excess death percentage. Herein unintentional injury and chronic lower respiratory disease have on average the highest excess death percentage.

cod_data %>%
  na.omit %>%
  group_by(cause_of_death) %>%
  mutate(mean_percent_excess_death = mean(percent_excess_death)) %>%
  ungroup(cause_of_death) %>%
  mutate(cause_of_death = fct_reorder(cause_of_death, mean_percent_excess_death)) %>%
  ggplot(aes(x = cause_of_death, y = percent_excess_death, fill = cause_of_death)) +
  geom_boxplot() +
  facet_grid(~locality) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))  +
  theme(legend.position = "bottom")

cod_lm_yz1 = lm(percent_excess_death ~cause_of_death, data = cod_data)
broom::tidy(cod_lm_yz1)

##                                              term estimate  std.error
## 1                                     (Intercept) 21.60686 0.08399078
## 2 cause_of_deathChronic Lower Respiratory Disease 21.08880 0.12091802
## 3                     cause_of_deathHeart Disease 10.89773 0.11879905
## 4                            cause_of_deathStroke 12.92196 0.12046778
## 5              cause_of_deathUnintentional Injury 25.77100 0.11875825
##   statistic p.value
## 1 257.25277       0
## 2 174.40579       0
## 3  91.73246       0
## 4 107.26483       0
## 5 217.00389       0

#PLOTLY

map_cod_data = cod_data %>%
  filter(locality == "Metropolitan") %>%
  select(state, locality, percent_excess_death) %>% 
  group_by(state) %>% 
  summarise(mean_ped = mean(percent_excess_death)) %>% 
  dplyr::filter(!(state == "District of\nColumbia"))
  

map = as.tibble(fifty_states) %>%
  group_by(id) %>% 
  summarize(clong = mean(long), clat = mean(lat)) %>% 
  filter(!(id == "district of columbia"))

df <- cbind(map, state.abb, state.center, rate = unique(map_cod_data$mean_ped))

ggplot(df, aes(map_id = id)) + 
    geom_map(aes(fill = rate), map = fifty_states) + 
    expand_limits(x = fifty_states$long, y = fifty_states$lat) + 
    labs(x = "", y = "") +
    theme(panel.background = element_blank(), 
          axis.text.x = element_blank(), 
          axis.text.y = element_blank(), 
          axis.ticks = element_blank()) + 
    geom_text(aes(x = clong, y = clat, label = state.abb)) +
  scale_fill_gradient(low="gold", high="red")

region_cod_data = cod_data %>%
  select(state, locality, hhs_region, percent_excess_death) %>% 
  group_by(state,locality, hhs_region) %>% 
  summarise(mean_ped = mean(percent_excess_death)) %>% 
  dplyr::filter(!(state == "District of\nColumbia")) %>% 
  mutate(hhs_region = as.character(hhs_region))


region_lm = lm(region_cod_data$mean_ped~region_cod_data$hhs_region)
rl1 = broom::tidy(region_lm)
kable(rl1)

region_locality_lm =lm(region_cod_data$mean_ped ~ region_cod_data$hhs_region+region_cod_data$locality)
rl2 = broom::tidy(region_locality_lm)
kable(rl2)